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Abstract 

Viscoelastic  materials  with  fading  memory,  e.g.,  polymers,  suspensions,  and  emulsions, 
exhibit  behavior  that  is  intermediate  between  the  nonlinear  hyperbolic  response  of  purely 
elastic  materials  and  the  strongly  diffusive,  parabolic  response  of  viscous  fluids.  Many 
popular  numerical  methods  used  in  the  computation  of  steady  viscoelastic  flows  fail  in 
important  flow  regimes,  and  thus  do  not  capture  significant  non-Newtonian  phenomena. 

A  key  to  satisfactory  explanation  of  these  phenomena  is  the  study  of  the  full  dynamics  of 
the  flow. 

This  paper  studies  the  dynamics  of  shear  flow,  presenting  a  description  of  non-Newtonian 
phenomena  caused  by  a  non-monotone  relation  between  the  steady  shear  stress  and  shear 
strain-rate.  Analytical  results  for  such  phenomena  are  surveyed,  and  three  distinct  nu¬ 
merical  methods  are  developed  to  accurately  compute  the  dynamics.  The  computations 
reproduce  experimental  measurements  of  non-Newtonian  “spurt”  in  shearing  flow  through 
a  slit  die.  They  also  predict  related  phenomena  (such  as  hysteresis  and  shape  memory); 
experiments  are  suggested  to  verify  these  predictions. 
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1.  Introduction 

Viscoelastic  materials  with  fading  memory,  e.g .,  polymers,  suspensions,  and  emulsions, 
exhibit  behavior  that  is  intermediate  between  the  nonlinear  hyperbolic  response  of  purely 
elastic  materials  and  the  strongly  diffusive,  parabolic  response  of  viscous  fluids.  Their 
properties  reflect  a  subtle  dissipative  mechanism  induced  by  the  fading  memory.  Several 
interesting  physical  phenomena,  which  are  important,  for  example,  in  polymer  processing, 
arise  in  shear  flows  of  viscoelastic  fluids.  Understanding  such  phenomena  has  proved 
to  be  of  significant  physical,  mathematical,  and  computational  interest.  We  have  found 
that  satisfactory  explanation  and  modeling  requires  the  study  of  the  full  dynamics  of  the 
equations  of  motion  and  the  constitutive  assumptions. 

One  striking  phenomenon  has  been  observed  by  Vinogradov,  et  a  1.  [18]  in  the  flow 
of  viscoelastic  fluids  (monodisperse  polyisoprenes)  through  capillaries.  They  found  that 
the  volumetric  flow  rate  increased  dramatically  at  a  critical  stress  that  was  independent 
of  molecular  weight.  This  phenomenon,  which  is  called  “spurt,”  had  been  overlooked 
or  dismissed  by  rheologists  because  no  plausible  mechanism  was  known  to  explain  it  in 
the  context  of  steady  flows.  Spurt  was  lumped  together  with  instabilities  such  as  “slip,” 
“apparent  slip,”  and  “melt  fracture,”  which  are  poorly  understood.  While  regarded  as 
anomalous,  these  instabilities  can  severely  disrupt  polymer  processes;  they  can  be  avoided 
in  practice  only  with  ad  hoc  engineering  expedients.  The  mechanisms  of  such  phenomena 
are  not  understood;  this  is  because  the  governing  equations  are  analytically  intractable,  and 
because  many  numerical  methods  for  steady  viscoelastic  fluid  flows  falter  in  this  regime, 
and  thus  cannot  model  the  spurt  phenomenon. 

Several  explanations  have  been  offered  for  the  spurt  phenomenon  [2,  4,  9,  12].  Their 
common  feature  is  that  the  shear  stress  in  steady  flow  does  not  vary  monotonically  with 
shear  strain  rate  (as  illustrated  in  Fig.  2,  below).  These  explanations  have  been  rejected 
by  many  rheologists  as  being  somehow  unphysical.  We  believe  that  this  criticism  is  un¬ 
founded  because  it  is  based  on  intuition  derived  from  generalized  Newtonian  models  of 
non-Newtonian  fluids. 

A  key  to  understanding  the  spurt  phenomenon  is  the  dynamical  behavior  of  the  con¬ 
stitutive  relations  as  well  as  the  equations  of  motion.  While  there  is  a  great  variety  of 
constitutive  models  for  viscoelastic  fluids,  the  dynamical  behavior  for  many  is  difficult  to 
analyze  or  compute.  In  this  paper,  we  model  the  spurt  phenomenon  using  the  Johnson- 
Segalman  model  [7]  of  a  non-Newtonian  fluid.  This  constitutive  relation  correctly  models 
the  spurt  phenomenon  and  yet  is  sufficiently  simple  to  be  understood  through  a  combina¬ 
tion  of  analysis,  asymptotics,  and  numerical  simulation.  A  non-monotone  stress-strain-rate 
relation  of  the  kind  that  causes  the  spurt  phenomenon  arises  when  the  fluid  behavior  is 
characterized  by  multiple  relaxation  times.  Interpretation  of  small-amplitude  oscillatory 
shear  data  in  Ref.  [18]  indicates  that  the  relaxation  times  are  widely  spaced.  Formal 
asymptotic  analysis  [10]  of  the  dynamics  shows  that  the  effects  of  the  smallest  relaxation 
time  are  mimicked  by  a  Newtonian  viscosity  term.  For  simplicity,  then,  we  study  the 
Johnson-Segalman  model  with  a  single  relaxation  time  and  added  Newtonian  viscosity. 

We  study  idealized  shearing  flow  through  a  narrow  slit  die.  Assuming  that  the  driving 
pressure  is  transmitted  instantaneously,  the  three-dimensional  flow  may  be  approximated 
by  a  one-dimensional  problem.  Our  analytical  and  numerical  results  show  that  flow  in  a  slit 
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die  reflects  the  essential  features  observed  for  capillaries.  We  believe  that  this  is  because 
the  spurt  phenomenon  depends  solely  on  material  properties  and  the  smallest  physical 
dimension  of  the  problem. 

The  outline  of  this  paper  is  as  follows.  In  §2,  we  describe  the  modeling  of  shearing 
flow  of  a  viscoelastic  fluid  using  the  Johnson-Segalman  constitutive  relation.  Mathematical 
results  for  this  model,  as  well  as  for  related  models  that  capture  some  key  features,  axe 
surveyed  in  §3.  These  include:  evolutionarity  of  the  system;  existence  and  regularity  of 
solutions;  formation  of  discontinuities;  asymptotic  behavior  for  large  time;  stability  of 
steady  solutions;  structure  of  discontinuous  solutions;  and  dynamics  of  a  related  system 
of  ordinary  differential  equations.  In  §4,  we  develop  three  distinct  numerical  methods  for 
solving  the  flow  equations,  accounting  for  the  mathematical  structure  of  the  model.  That 
these  methods  reproduce  physical  phenomena  is  demonstrated  in  §5,  where  we  compare 
numerical  calculations  with  experimental  data  for  the  spurt  phenomenon.  Based  on  these 
results,  we  propose  in  §6  some  rheological  experiments  to  confirm  the  predictions  of  the 
model.  Finally,  in  §7,  we  discuss  our  conclusions. 

2.  Mathematical  Formulation 

2a.  Shear  Flow  of  a  Johnson-Segalman  Fluid 

The  motion  of  a  fluid  under  incompressible  and  isothermal  conditions  is  governed  by 
the  balance  of  linear  momentum 

=  VS.  (2.1) 

Here,  p  is  the  fluid  density,  v  is  the  particle  velocity,  and  S  is  the  stress  tensor.  The 
response  characteristics  of  the  fluid  are  embodied  in  the  constitutive  relation  for  the  stress. 
For  viscoelastic  fluids  with  fading  memory,  these  relations  specify  the  stress  as  a  functional 
of  the  deformation  history  of  the  fluid.  Many  sophisticated  constitutive  models  have  been 
devised;  see  Ref.  [1]  for  a  survey.  In  the  present  work,  we  focus  on  the  Johnson-Segalman 
model  [7]  as  a  prototype  for  general  constitutive  models.  This  model  accounts  for  non- 
affine  deformation  of  Gaussian  networks  by  introducing  a  slip  parameter  a,  —  1  <  a  <  1, 
leading  to  a  nonlinear  generalization  of  the  classical  Maxwell  model. 

To  specify  this  constitutive  relation,  we  decompose  the  stress  as 

S  =  -pi  +  2pD  4-  E  .  (2.2) 

In  this  equation,  p  is  an  isotropic  pressure  (which  is  determined  from  the  incompressibility 
constraint),  rj  is  the  coefficient  of  Newtonian  viscosity,  and  E  is  the  non-Newtonian  extra 
stress.  Also,  we  let  D  :=  \  [Vv  +  (Vv)T]  and  fi  :=  \  [Vv  -  (Vv)r]  be  the  symmetric 
and  antisymmetric  parts  of  the  velocity  gradient  Vv,  which  has  components  (Vv)'j  := 
dv'  jdxi .  The  extra  stress  is  specified  by  the  differential  constitutive  law 


dv 

dt 


4.  v  •  Vv 


E  =  2pB  -  AS  ,  (2.3) 
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where 

£  :=  ^  +  v  •  VS  +  E[fl  -  dD]  +  [O  -  aD]TE  (2.4) 

is  the  objective  time  derivative  of  E  with  parameter  a.  The  parameter  p  is  an  elastic  shear 
modulus,  and  A  is  a  relaxation  rate. 

Constitutive  relations  such  as  Eq.  (2.3)  exhibit  a  mixture  of  elastic  and  viscous  be¬ 
havior.  This  may  be  seen  heuristically  as  follows.  In  the  long  relaxation-time  limit,  A  — ►  0, 

Eq.  (2.3)  shows  that  an  objective  time  derivative  of  E  is  proportional  to  the  deformation 

* 

rate:  E  ~  2pD.  This  is  characteristic  of  elastic  behavior,  and  leads  to  the  interpretation 
of  /x  as  a  shear  modulus.  By  contrast,  when  A with  nj\  fixed,  E  2(^/A)D; 
thus  the  model  displays  viscous  behavior  with  p/A  being  the  Newtonian  shear  viscosity 
coefficient. 

Essential  properties  of  the  constitutive  relation  are  exhibited  in  simple  planar  shear 
flow.  With  the  flow  aligned  along  the  y-axis  (see  Fig.  1),  the  flow  variables  are  independent 
of  y.  Therefore  the  velocity  field  is  v  =  (0,  v(x ,  f)),  and  the  balance  of  mass  is  automatically 
satisfied.  Furthermore,  the  components  of  the  extra  stress  tensor  E  may  be  written  E11  = 
7(1,  t),  Eiy  =  E91  =  a(x,  t),  and  Ey!/  =  r(r,  f),  while  the  pressure  takes  the  form  p  = 
po(ar,  t )  —  f  being  the  pressure  gradient  driving  the  flow.  In  these  terms,  Eqs.  (2.3) 

become 


7t  +  (1  ~  a)<™x  =  -A7  , 

(2.5  a) 

at  -  [|(1  +  a) 7  -  |(1  -  a)r  +  fj]  vx  =  -A a 

,  (2.56) 

rt  —  (1  +  a)avx  =  -Ar  . 

(2.5c) 

variables  Z  :=  j(l-t-c)7—  |(1  —a)r  and  W  := 

—  j(l  +0)7—  A(l.  —  a)r, 

at  -{Z  +  n)vx  —  -A a  , 

(2.6a) 

Zt  +  (1  —  a2)avx  =  —A Z  , 

(2.66) 

Wt  =  -A W  . 

(2.6c) 

Because  W  must  remain  finite  as  t  — ►  —00,  W  =  0,  and  the  last  equation  may  be  omitted. 
As  a  result,  Z  =  —  ~(1  —  a2)(r  —  7),  where  Eyy  —  E1*  =  r  —  7  is  the  principal  normal  stress 
difference. 

Combining  the  constitutive  law  (2.6)  with  the  balance  of  linear  momentum  (2.1),  we 
are  led  to  the  system  of  equations 


pvt  -  <7i  =  T]vxx  +  f  ,  (2.7 a) 

at  -  (Z  +  n)vx  =  -A a  ,  (2.76) 

Zt  +  (1  —  a2)avx  =  —A Z  .  (2.7c) 

In  this  paper,  we  study  shear  flow  between  two  parallel  plates,  located  at  x  =  ±h/2.  By 


symmetry,  we  need  only  consider  the  flow  on  the  interval  [—h/ 2,0].  The  no-slip  condition 
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Fig.  1:  Shear  flow  through  a  slit-die. 

at  the  plate  implies  the  boundary  condition  v(—h/ 2,  t)  =  0,  while  symmetry  imposes  that 
vx(0 ,t)  =  0.  We  also  prescribe  initial  values  for  v,  cr,  and  Z,  which  must  be  compatible 
with  the  boundary  conditions.  To  conform  with  the  symmetry,  we  require  that  a( 0, 0)  =  0; 
then,  according  to  Eq.  (2.7b),  cr(0,  t)  =  0  for  all  time. 

To  eliminate  unnecessary  parameters,  we  scale  distance  by  h,  time  by  A-1,  and  the 
stresses  a  and  Z  by  p.  Furthermore,  if  we  replace  cr,  v,  and  /  by  a  :=  (1  —  a2)1/2^, 
v  :=  (1  —  a2y/2v,  and  /  :=  (1  —  a2)1/2/,  respectively,  then  the  parameter  a  disappears  from 
Eqs.  (2.7).  Since  no  confusion  will  arise,  we  omit  the  caret.  The  dimensionless  parameters 
are  a  :=  ph2 A2/^i  and  e  :=  tj X/p.  Consequently,  we  study  the  initial-boundary-value 
problem  for  the  system 

otvt  —  ax  =  evxx  4-  /  , 

-  (Z  +  l)vx  =  -cr  ,  (JS) 

Zt  +  crvx  =  —Z  , 

on  the  interval  [-1/2,0],  with  boundary  conditions 

t>(  — 1/2,  t)  =  0  and  vx(0,t)  =  0  {BC) 
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and  initial  conditions 


u(x,0)  =  v0(x )  ,  <r(x,0)  =  <70(x)  ,  and  Z(x,0 )  =  Z0(x )  ,  (IC) 

where  t>o(— 1/2)  =  0,  Vq(0 )  =  0  and  Co(0)  =  0.  For  later  purposes,  notice  that  the 
momentum  equation  may  be  written  a  vt  —  Tx  =  f ,  where  T  :=  a  +  evx  denotes  the  total 
shear  stress. 


2b.  Steady  Shear  Flow 

The  steady-state  solutions  of  (JS),  when  the  forcing  term  /  is  a  constant  /,  play  an 
important  role  in  our  discussion.  Such  a  solution,  denoted  by  u,  a,  and  Z,  is  given  as 
follows.  The  stress  components  o’  and  Z  are  related  to  the  velocity  gradient  vx  (which,  in 
dimensionless  units,  is  the  Deborah  number)  through 


a  = 


Vx 


(2.8) 


and 


Z  +  l  = 


1 

l  +  vl  ' 


(2.9) 


Therefore  the  total  steady  shear  stress  T  :=  cr  +  evx  is  given  by  T  =  Tstea.dy(vx),  where 


Tsteady(ux)  :==  d-  CU 

1  +vi 


(2.10) 


When  €  <  1/8,  this  relation  between  the  steady  shear  stress  and  strain  rate  in  not  mono¬ 
tone,  as  illustrated  in  Fig.  2.  In  this  manner,  the  “stress  loop”  of  Fig.  2  arises  automatically 
in  the  Johnson- Segalman  model. 

The  momentum  equation,  together  with  the  boundary  condition  at  the  centerline, 
implies  that  the  total  steady  shear  stress  satisfies  T  =  —fx  for  x  6  [—  j,0].  Therefore  the 
velocity  gradient  may  be  determined  as  a  function  of  x  by  solving 

Tsteady(^x)  =  —  fx  .  (2.11) 


The  steady  velocity  profile,  shown  in  Fig.  3,  is  obtained  by  integrating  vx  and  using  the 
boundary  condition  at  the  wall.  However,  because  the  function  Tsteady  is  not  monotone, 
there  may  be  up  to  three  distinct  values  of  vx  that  satisfy  Eq.  (2.11)  for  any  given  x. 
Consequently,  vx  may  suffer  jump  discontinuities,  resulting  in  kinks  in  the  velocity  profile 
(as  at  the  point  x,  in  Fig.  3).  Indeed,  a  steady  solution  must  contain  such  a  jump  if  the 
total  stress  Twaii  =  J/2a.t  the  wall  exceeds  the  total  stress  at  the  local  maximum  in  Fig.  2; 
for  later  convenience,  we  denote  this  critical  stress  by  Tcrit  =  2- 
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^wall 

Tent 

-/z* 


O 


Fig.  2:  Total  steady  shear  stress  T  vs.  shear  strain  rate  vx  for 
steady  flow. 

2c.  Comparison  with  the  Generalized  Newtonian  Model 

Traditionally,  a  non-monotone  relation  between  stress  and  strain  rate  is  regarded  as 
a  defect  of  the  constitutive  law.  This  conclusion  is  based  on  intuition  appropriate  for 
generalized  Newtonian  models  of  non-Newtonian  fluids.  Shear  flow  for  such  a  fluid  is 
governed  by  the  single  equation 


pvt  -  [tj(v *)i7r]„  =  /  ,  (2.12) 

corresponding  to  having  a  viscosity  coefficient  77  that  depends  on  strain-rate.  In  a  flow 
regime  where  t]{vx)vz  decreases  with  strain  rate  vx,  however,  Eq.  (2.12)  has  the  character 
of  a  backward  heat  equation,  which  suffers  from  the  Hadamard  instability.  Therefore  for 
generalized  Newtonian  fluids,  tj{vx)vx  must  increase  with  vx  in  a  physically  stable  steady 
solution. 

The  system  (JS)  has  the  same  steady  solutions  as  a  generalized  Newtonian  model  with 
T}{vx)vx  ~  Tsteady(vx),  so  one  might  think  that  it  exhibits  the  same  instability  in  regions 
where  Tsteady  decreases.  This  conclusion  is  not  warranted,  however,  because  the  system 
(JS)  maintains  its  evolutionary  character  when  e  >  0  (see  §3a). 


-7- 


Non-Newtonian  Shear  Flow 


Malkus,  Nohel,  Plohr 


Fig.  3:  Velocity  profile  for  steady  flow. 

3.  Mathematical  Results 

Several  mathematical  results  are  known  for  the  system  ( JS);  we  refer  to  Refs.  [6,  17, 

4,  16,  3,  14,  10]  for  further  discussions  and  additional  references. 

3a.  Evolutionarity,  Existence  of  Solutions,  and  Formation  of  Discontinuities 
When  the  viscosity  parameter  e  =  0,  the  quasi-linear  system  (JS)  is  strictly  hyperbolic 
provided  that  Z  4-  1  >  0.  In  this  case,  the  wave  speeds  are  db  [(Z  4-  l)/a]1//2  and  zero.  If, 
on  the  other  hand,  Z  +  1  becomes  negative,  then  (JS),  with  e  =  0,  undergoes  a  change  of 
type  and  loses  its  evolutionary  character.  Joseph,  Renardy,  and  Saut  [6]  have  associated 
this  change  of  type  with  certain  fluid  instabilities. 

Suppose  that  e  =  0  and  /  =  0,  and  assume  that  the  initial  data  are  smooth  and  lie  in 
the  hyperbolic  region.  If  the  data  have  sufficiently  small  variation,  then  a  unique  classical 
solution  of  (JS),  (IC),  (BC)  exists  globally  in  time;  moreover,  the  solution  decays  to  zero 
as  t  — ►  oo.  This  can  be  proved  using  the  energy  methods  discussed  in  Ref.  [17].  On  the 
other  hand,  if  the  data  have  sufficiently  large  variation,  then  the  classical  solution  blows 
up  within  finite  time;  i.e.,  |uz|,  jcrx |,  and  |z*|  approach  infinity  as  t  approaches  a  finite 
critical  time.  This  is  proved  in  Ref.  [17]  using  the  method  of  characteristics.  Thus  the 
fading  memory  acts  as  a  weak  dissipative  mechanism:  the  source  terms  in  the  equations 
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serve  to  counteract  the  formation  of  singularities  from  sufficiently  smooth  data.  When 
discontinuities  do  form,  system  (JS)  is  no  longer  valid  because  the  products  of  distributions 
Zvx  and  cvx  are  ill-defined.  (See  the  discussion  in  §3b.) 

If  e  >  0,  the  system  (JS)  is  evolutionary,  but  it  cannot  be  classified  according  to 
type.  Recently,  Guillope  and  Saut  [3]  established  the  global  existence  of  solutions  of  (JS) 
for  planar  Couette  and  Poiseuille  flow  with  data  of  arbitrary  size.  They  also  studied  the 
asymptotic  (Lyapunov)  stability  of  steady  states  in  the  Couette  case. 

3b.  Conservation  Laws 

It  is  important  to  observe  that  (JS)  is  not  in  conservation  form.  The  evolution  of  a 
Johnson-Segalman  fluid  is,  in  fact,  governed  by  physical  conservation  laws  [7].  A  conser¬ 
vative  formulation  of  (JS)  must  be  used  when  <r,  u,  and  Z  are  discontinuous. 

Following  Plohr  [16],  we  introduce  the  “elastic  part”  t  of  the  shear  strain  and  the 
“entropy”  variable  z  through  the  relations 

a  ~  z  sinr  , 

Z  +  1  :=  z  cos  t  . 

Then  system  (JS)  is  transformed  into  the  equivalent  system 

rt  —  vx  —  —z~1  sin r  , 
avt  -  [<r(r,z)  +  evx}x  =  /  , 

zt  =  — (z  —  cost)  , 

which  is  in  conservative  (i.e.,  divergence)  form.  Furthermore,  if  the  internal  energy  £  is 
defined  by 

a  £  :=  1  —  z  cos  r  ,  (3.2) 

the  energy  is  dissipated  according  to  the  equation 

a  [iu2  +  £(t,  z)\  t  -  {[<r(r,  z)  +  e  vx]  v}t  =  vf  -  a  £ (r,  z)  -  e  ( vx  f  .  (3.3) 

The  conservative  formulation  (C)  of  (JS)  is  used  in  the  numerical  methods  discussed  in 
§4c. 


(C) 


(3-la) 

(3.16) 
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3c.  Model  Problems 

More  detailed  analytical  results  are  obtained  by  simplifying  the  system  (JS).  A  model 
system  that  incorporates  several  qualitative  features  of  (JS)  is  obtained  by  freezing  Z  at 
its  equilibrium  value:  Z  +  1  =  1/(1  +  v\ ).  Defining  g(vx )  :=  vx/(l  +  v%),  system  (JS) 
becomes 

avt  -  <rx  =  evxx  +  f  , 

r  s  (M) 

~  g{vx)  =  -cr  . 

More  generally,  g  may  be  any  smooth,  odd  function.  The  boundary  and  initial  conditions 
for  v  and  a  are  the  same  as  in  (BC)  and  (IC).  We  assume  that  e  >  0  and  that  /  is  the 
constant  /.  The  function  g  is  related  to  the  steady  stress-strain-rate  relation  through 
Tsteady(vx)  =  g(v x)  +  tvx.  A  steady  solution  of  (M)  satisfies  W  =  g(vx )  and  Tste ady(v*)  = 
— fx,  just  as  for  the  system  (JS). 

Nohel,  Pego,  and  Tzavaras  [14]  have  shown  that  the  global  classical  solution  v,  a  of 
(M),  (BC),  (IC)  has  the  following  properties. 

(i)  With  S  :=  cr  +  evx  -(-  /x,  S(x,t)  — ►  0  as  t  — >  oo,  uniformly  for  x  £  [—1/2,0]. 

(ii)  There  exists  a  steady  state  u,  or  such  that  for  each  x  £  [—1/2,0],  v(x,t )  — ►  v(x), 
vx(x,t)  —*  vx{x),  and  cr(x,t)  — >  c(x)  as  t  — ►  oo.  We  emphasize  that  the  steady  velocity 
gradient  vx  and  stress  <F  may  be  discontinuous  (as  in  Fig.  2). 

(iii)  Let  u,  cr  be  a  steady  state  such  that 

^'teadyO7*)  =  9'(vx)  4-  e  >  const.  >  0  .  (3.4) 

(Referring  to  Fig.  2,  inequality  (3.4)  precludes  top  and  bottom  jumping  and  excludes  the 
region  where  Tsteady(ui)  decreases.)  Consider  a  union  U  of  small  subintervals  of  —  |  < 
x  <  0  that  are  centered  at  points  where  vx  and  a  are  discontinuous.  Let  smooth  initial 
data  be  chosen  such  that  |S(x,0)|  is  sufficiently  small  except  in  U.  Then  the  solution  of 
(M)  converges  to  the  steady  state  v,  a  on  the  complement  of  U.  Moreover,  the  measure  of 
U  can  be  made  arbitrarily  small  by  choosing  |S(x,0)|  small  enough.  In  this  sense,  steady 
states  are  stable  (even  if  vx  and  <?  are  discontinuous). 

The  numerical  results  discussed  in  §5  suggest  that  similar  results  hold  for  the  system  (JS). 
Proofs  for  (JS)  are  under  investigation. 

The  model  problem  (M)  was  studied  also  by  Hunter  and  Slemrod  [4],  In  their  con¬ 
struction  of  the  model,  the  steady-state  relation  W  =  g(vx )  between  the  stress  and  strain 
rate  is  chosen  to  be  gtis(vx)  •'=  chs{vx)  —  evx,  where  the  graph  of  the  function  chs  resem¬ 
bles  Fig.  2.  In  contrast  with  ojs  '•=  Tsteady,  cths  is  independent  of  e.  Hunter  and  Slemrod 
base  their  analysis  on  the  conservation  laws 


Wt  ~  ux  —  0  ,  (3.5a) 

aut  -  chs{w)z  =  euxx  —  au  (3.5 b) 

for  the  acceleration  u  =  vt  and  the  strain  rate  w  =  vx.  Therefore  jumps  in  the  strain  rate 
vx  are  seen  to  correspond  to  steady  shock  waves  for  the  system  (3.5)  with  e  =  0.  Based 
on  a  local  dynamical  analysis  of  shock  structure  for  small  e,  the  centerline  velocity  is 
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shown  to  exhibit  hysteresis  under  quasi-static  cycling  of  the  pressure  gradient.  (This  same 
behavior  is  observed  in  the  numerical  simulation  of  the  system  (JS);  see  §5.)  We  emphasize, 
however,  that  this  analysis  cannot  be  applied  to  the  model  problem  (M)  as  derived  from 
the  Johnson-Segalman  system  (JS)  because  the  function  gjs(vx)  :=  ux/(  1  +  t>2)  decays  to 
zero  at  high  strain-rate. 

3d.  Phase-Plane  Analysis 

Quite  detailed  information  about  the  structure  of  solutions  of  (JS)  can  be  garnered  by 
studying  a  system  of  ordinary  differential  equations  that  approximates  it.  Motivation  for 
this  approximation  comes  from  the  following  observation:  in  experiments  of  Vinogradov, 
et  a  1.  [18],  a  =  ph2\2/p  is  of  the  order  10-12;  thus  the  term  avt  in  the  momentum 
equation  of  (JS)  is  negligible  even  when  is  moderately  large.  We  are  led  to  study  the 
approximation  to  (JS)  obtained  when  a  =  0.  As  we  now  outline,  the  behavior  of  solutions 
of  the  resulting  dynamical  system  offers  an  explanation  for  several  features  of  the  flows 
calculated  for  the  full  system  (see  §5);  in  fact,  these  calculations  prompted  the  following 
analysis.  A  detailed  exposition  of  these  results  is  to  be  found  in  Ref.  [10]. 

When  a  =  0,  the  momentum  equation  may  be  integrated,  just  as  in  the  case  of  steady 
flows,  to  show  that  the  total  shear  stress  T  :=  a  -(-  evx  coincides  with  the  steady  value 
T(x)  =  —fx.  The  remain;ug  equations  of  (JS)  become,  for  each  fixed  x,  the  autonomous 
planar  system  of  ordinary  differential  equations 

crt  =  (Z  +  1)  ~  a  >  (3.6a) 

-  Z  .  (3.6b) 

We  emphasize  that  a  different  dynamical  system  is  obtained  at  each  point  in  the  channel. 
These  dynamical  systems  can  be  analyzed  completely  by  a  phase-plane  analysis  [10],  with 
the  following  result. 

(i)  A  critical  point  for  the  system  (3.6)  defines  a  steady-state  solution  of  (JS);  such  a 
solution  corresponds  to  a  point  on  the  steady  total-stress  curve  (see  Fig.  2)  at  which  the 
total  stress  is  T(x).  Consequently,  there  is  a  single  critical  point  when  T  is  sufficiently 
small  or  sufficiently  large,  while  there  are  three  critical  points  when  T  is  intermediate  in 
value.  A  critical  point,  such  as  A,  that  lies  between  the  origin  O  and  the  local  maximum 
M  of  the  total-stress  curve  is  an  attracting  node,  which  we  call  a  classical  attractor.  A 
point,  such  as  C,  that  lies  past  the  local  minimum  m  is  either  an  attracting  node  or  an 
attracting  focus,  called  a  spurt  attractor.  A  critical  point  B  lying  between  M  and  m  is 
a  saddle  point.  The  phase  plane  for  an  example  with  three  critical  points  is  shown  in 
Fig.  4,  where  the  invariant  manifolds  for  the  classical  attractor  A  and  the  saddle  point  B 
have  been  drawn.  By  constructing  an  invariant  region  for  system  (3.6)  and  analyzing  its 
critical  points  at  infinity,  the  qualitative  structure  of  orbits  for  the  dynamical  system  can 
be  determined  completely. 

(ii)  Suppose  that  the  initial  conditions  for  the  flow  are  at  the  origin,  <7  =  0  and  Z  =  0,  and 
suppose  that  the  forcing  /  is  supercritical.  For  points  x  near  the  centerline,  where  T(x) 
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Fig.  4:  Phase  plane  when  T  is  subcritical.  The  critical  points  A, 

B ,  and  C  are,  respectively,  an  attracting  node,  a  saddle  point,  and 
a  stable  focus. 

lies  below  Tcriti  the  origin  lies  in  the  basin  of  attraction  of  the  classical  attractor.  This  is 
illustrated  in  Fig.  4,  where  the  integral  curve  starting  at  the  origin  O  ends  at  A.  When 
T(x )  exceeds  Tcriti  the  only  critical  point  is  the  spurt  attractor.  The  orbit  through  the 
origin  O  in  this  case  is  shown  in  Fig.  5.  Consequently  the  flow  is  predicted  to  approach  a 
steady  spurt  solution  in  which  the  jump  in  strain  rate  occurs  at  the  maximum  stress  (“top 
jumping”),  with  the  kink  in  the  velocity  profile  located  as  close  as  possible  to  the  wall. 
Similar  arguments  explain  the  hysteresis  effects  that  occur  upon  unloading. 

(iii)  More  quantitative  information  is  obtained  when  e  is  small.  Referring  to  Fig.  2,  the 
total  stress  Tcr it  at  the  the  local  maximum  M  is  1/2  0(e),  while  the  local  minimum  m 

corresponds  to  a  tc  .  u  stress  of  2>/e  [1  -j-  0(e)].  Furthermore,  cr  =  T  +  O(e)  at  a  classical 
critical  point,  while  a  —  0(e)  for  a  spurt  attractor.  Consider  a  point  along  the  channel  for 
which  T(x)  >  Tcriti  so  that  the  only  critical  point  is  the  spurt  attractor,  and  suppose  that 
that  T  <  1.  Then  the  evolution  of  the  system  exhibits  three  distinct  phases,  as  indicated 
in  Fig.  5:  an  initial  “Newtonian”  phase  (O  to  N)\  an  intermediate  “latency”  phase  ( N  to 
5);  and  a  final  “spurt”  phase  (5  to  O). 

The  Newtonian  phase  occurs  on  a  time  scale  of  order  e,  during  which  the  system 
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Fig.  5:  Phase  plane  when  T  is  supercritical.  The  point  C  is  the 
spurt  attractor;  point  L  is  located  at  Z  =  —  1  -f  T. 

approximately  follows  an  arc  of  a  circle  centered  at  a  —  0  and  Z  =  —  1.  Thus  Z  approaches 

^Newton  =  (1  -  T*)?  -1  (3.7) 

as  a  rises  to  the  value  T.  (If  T  >  1,  a  never  attains  the  value  T .)  The  latency  phase 
is  characterized  by  having  a  =  T  +  0(e),  so  that  a  is  nearly  constant  and  Z  evolves 
approximately  according  to  the  differential  equation 

Z'  =  -2TT-Z'  (3'8) 

Therefore  the  shear  stress  and  velocity  profiles  closely  resemble  those  for  a  steady  solution 
with  no  spurt,  but  the  solution  is  not  truly  steady  because  the  normal  stress  difference  Z 
still  changes.  Integrating  Eq.  (3.S)  from  Z  =  ^Newton  to  Z  =  -1  determines  the  latency 
period.  This  period  becomes  indefinitely  long  when  the  forcing  decreases  to  its  critical 
value;  thus  the  persistence  of  the  near-steady  solution  with  no  spurt  can  be  very  dramatic 
(see  §5).  The  solution  remains  longest  near  point  L  where  Z  =  —  1  +  T.  This  point  may  be 
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regarded  as  the  remnant  of  the  classical  attractor  A  and  the  saddle  point  B.  Eventually 
the  solution  enters  the  spurt  phase  and  tends  to  the  spurt  attractor.  Because  this  critical 
point  is  a  focus,  the  stress  oscillates  between  the  shear  and  normal  components  while  it 
approaches  the  steady  state. 


4.  Numerical  Methods 

To  study  the  dynamics  of  system  (JS),  we  developed  several  different  numerical  meth¬ 
ods;  each  has  its  advantages  for  certain  ranges  of  physical  parameters.  Calculations  with 
these  methods  produce  similar  qualitative  and  quantitative  results. 

4a.  Solid  Mechanics  Formulation 

In  the  solid  mechanics  formulation,  the  system  (JS)  is  regarded  as  governing  the 
extensional  motion  of  an  elastic-plastic  bar.  The  first  equation  is  momentum  balance,  in 
which  the  parabolic  term  adds  viscous  “stiffness  damping.”  The  remaining  equations  are 
incremental  constitutive  relations  for  the  stress.  The  stiffness  of  the  material  is  reflected 
in  the  wave  speed  [( Z  +  1  )/<*]1,/2.  We  have  observed  that  the  wave  speed  is  diminished 
under  loading,  so  that  the  material  exhibits  plastic  softening.  (See  also  Ref.  [16]  for  an 
interpretation  of  (JS)  as  governing  a  viscoplastic  material.) 

We  have  solved  the  system  (JS)  numerically  using  a  method  motivated  by  solid  me¬ 
chanics.  This  method  is  based  on  a  Galerkin  weak  form  of  the  momentum  equation,  with 
test  functions  <j>  £  $  and  trial  functions  v(-,t)  £  $  for  each  t.  We  use  continuous  piecewise 
linear  functions  for  and  we  represent  a  and  Z  as  piecewise  constant.  The  semidiscrete 
Galerkin  equation  (integrated  by  parts)  is  evaluated  at  tn+i  =(n  +  l)At,  yielding 

f° 

/  {a<j>  (ut)n+1  -I-  <f>x<?n+i  +e<i>x  Mn+1  -  <f>fn+i}  dx  =  0  ,  (4.1) 

J —1/2 

where  (u£)n+1  means  v£(x,tn+i),  etc.  The  Galerkin  equation  is  solved  for  vn+i  by  ad¬ 
vancing  the  viscoelastic  contribution  to  the  stress  using  a  semi-implicit  treatment  of  the 
constitutive  equation: 


crn+i  =  (1  -  At)an  +  AtZn(vz)n  +  At  (vr)n+1  •  (4.2) 

The  time  derivative  is  discretized  with  a  trapezoidal  approximation,  so  that 

«n+l  =  vn  +  At  {(1  -  7)  (v«)n  +  7  (t>0n+ 1  }  (4-3) 

for  some  parameter  7  >  0.  The  combination  of  Eqs.  (4.1)-(4.3)  is  a  method  in  which  the 
damping  term,  which  has  effective  viscosity  e-f  At,  is  treated  implicitly,  while  the  nonlinear 
term  involving  Zvx  is  treated  explicitly.  The  matrix  formulation  is  described  in  Ref.  [8];  as 
with  linear  problems,  the  system  matrix  needs  refactorization  only  if  At  is  changed.  After 
solution  of  the  combined  equations  for  u„+i,  o  is  corrected  and  Z  is  advanced  according 

<7n+i  =  (1  —  At)  crn  4-  At  (Z„  +  1)  (vr)n+1  ^  ^ 

Zn+l  =  (1  —  At)  Zn  —  At(Tn-(-i  (Vz)n+1 
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Other  variants  of  (4.4)  that  treat  the  constitutive  equations  with  a  higher  degree  of  im¬ 
plicitness  are  obvious.  However,  as  long  as  they  are  combined  with  (4.2),  which  limits  the 
scheme  to  first  order  accuracy  in  time,  there  is  no  reason  to  expect  such  variants  to  be 
more  accurate.  They  may  improve  stability,  but  the  algorithm,  as  given,  is  very  stable. 

The  stability  of  this  method  has  been  analyzed  for  the  system  (JS)  with  Z  frozen  [11]: 
the  method  is  stable  provided  that  Z  + 1  >  —  e  and  the  time  step  is  restricted  by  At  <  2 
(i.e.,  At  <  2/A,  dimensionally).  Of  course  if  Z  +  1  <  — e,  then  the  differential  equations 
themselves,  as  well  as  the  method,  are  linearly  unstable. 

4b.  Parabolic  Formulation 

Recall  that  the  total  stress  is  defined  to  be  T  =  a+ evx\  assume  that  e  >  0.  Introducing 
T  as  an  independent  variable,  the  system  (JS)  is  replaced  by 

Tt  =  ^Txx  +  (Z  + 1)  (^-~ :)  “  *  , 

**t  =  (Z  +  1)  f- — -'j  -  <y  , 

Zt  =  -*(lzz  )-z. 

The  boundary  conditions  are  Tx{— 1/2,  t)  =  —f  and  T(0,  t)  =  0.  The  velocity  profile  may 
be  reconstructed  by  integrating  (T  —  cr)/e. 

The  system  (4.5)  has  the  form  of  a  linear  heat  equation  forced  by  a  nonlinear  heat 
source  that  is  governed  by  two  auxiliary  ordinary  differential  equations.  To  solve  this 
system  numerically,  we  discretize  the  parabolic  term  in  Eq.  (4.5a)  implicitly  while  treating 
the  remaining  forcing  terms  explicitly.  Time  integration  is  performed  using  a  packaged 
ordinary  differential  equation  solver. 

We  remark  that  system  (4.5)  is  convenient  also  for  studying  existence  and  regularity 
of  solutions  of  (JS). 

4c.  Conservative  Formulation 

The  system  (JS)  is  equivalent  to  the  system  (C);  therefore  it  may  be  studied  from  the 
viewpoint  of  conservation  laws.  In  Ref.  [16],  we  have  determined  completely  the  structure 
of  scale-invariant  nonlinear  waves  for  (C)  when  e  =  0.  Such  a  wave  consists  of  a  sequence 
of  elementary  scale-invariant  waves,  either  centered  discontinuities  or  rarefaction  waves, 
connecting  constant  states  on  the  left  and  right.  Discontinuities  are  required  to  satisfy  Liu’s 
generalization  of  Oleinik’s  entropy  condition,  which  guarantees  that  energy  is  dissipated 
( cf.  Eq.  (3.3)).  This  admissibility  condition  is  equivalent  to  requiring  shock  waves  to  have 
viscous  profiles:  admissible  shock  waves  arise  as  limits  of  traveling-wave  solutions  of  (C) 
as  e  — +  0.  Our  \nalysis  follows  the  techniques  for  general  systems  of  conservation  laws 
discussed  in  Refs.  [13]  and  [5]. 

The  wave  structure  is  conveniently  depicted  with  a  wave  curve,  the  locus  of  states 
U  =  ( t,v,z )  on  the  right  for  a  fixed  state  Ui  =  ( tl,vl,zl )  on  the  left.  Only  z  may 
change  across  waves  with  zero  speed,  so  these  wave  curves  are  trivial.  On  the  other  hand, 
z  remains  constant  across  waves  corresponding  to  the  characteristic  families  with  speeds 


(4.5a) 

(4.56) 

(4.5c) 
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dbc,  so  we  suppress  the  z  component  of  U  when  drawing  wave  curves.  Figure  6  shows  the 
wave  curve  of  the  — c  (i.e.,  slow)  family  for  a  representative  initial  state  Ui;  through  points 
along  this  curve  are  drawn  the  wave  curves  of  the  +c  (i.e.,  fast)  family.  The  figure  was 
produced  using  a  computer  program  that  constructs  the  wave  curves  for  general  systems 
of  conservation  laws  [5]. 


Fig.  6:  The  wave  curve  diagram  for  a  typical  Ui.  Solid  curves  cor¬ 
respond  to  rarefaction  waves,  dashed  curves  to  shock  waves,  and 
dashed  curves  with  crossbars  to  composite  waves.  Darker  curves 
represent  the  — c  family,  while  lighter  curves  represent  the  -j-c  fam¬ 
ily. 

With  the  structure  of  scale-invariant  waves  known,  Riemann  initial-value  problems 
may  be  solved.  This  is  illustrated  in  Fig.  6  for  the  Riemann  problem  with  left  state  Ui 
and  right  state  Ur:  the  solution  contains  the  middle  state  Um  separating  a  slow  composite 
wave  on  the  left  from  a  fast  composite  wave  on  the  right.  We  have  written  a  computer 
program  that  solves  Riemann  problems,  and  have  incorporated  it  into  the  Glimm-Chorin 
random  choice  method.  This  method  solves  the  Cauchy  problem  without  introducing 
artificial  Newtonian  viscosity.  We  refer  to  Ref.  [15]  for  a  detailed  discussion. 
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5.  Numerical  Results 

In  this  section  we  describe  several  features  of  the  shear  flow  of  a  Johnson-Segalman 
fluid,  as  obtained  using  the  numerical  methods  of  §4. 

5a.  Effects  of  Newtonian  Viscosity 

As  our  first  numerical  experiment,  we  simulated  system  (C)  with  e  =  0  using  the 
random  choice  method.  Parameters  were  chosen  so  that  a  =  1.  The  flow  was  initially  in 
the  classical  steady  state  corresponding  to  the  critical  pressure  gradient  /crit  =  1;  then  the 
pressure  gradient  is  increased  abruptly  to  the  super-critical  value  1.2/crit. 

The  result  is  shown  in  Fig.  7.  The  fluid  velocity  v  is  plotted  vs.  position  x  at  successive 
time  intervals;  generally  the  velocity  increases  with  time.  During  the  early  stages  of  the  ex¬ 
periment,  the  flow  settled  into  a  quasi-steady  state.  This  latency  effect  is  especially  evident 
in  a  plot  of  the  centerline  velocity  as  a  function  of  time,  and  it  is  more  pronounced  when 
a  is  smaller.  Eventually,  however,  a  thin  layer  develops  at  the  plate  in  which  the  velocity 
rises  to  a  value  that  is  nearly  constant  across  the  channel.  For  practical  purposes,  the  fluid 
has  broken  free  from  the  plate  and  is  accelerating  uniformly  under  the  applied  pressure 
gradient;  thus  the  fluid  “slips.”  This  occurrence  might  be  related  to  the  phenomenon  of 
wall  slip,  which  has  been  associated  with  non-monotone  constitutive  relations  [2,  9]. 

It  is  worth  noting  that  the  random  choice  method  is  the  only  numerical  method 
that  can  calculate  super-critical  flows  with  e  =  0:  other  methods  must  be  stabilized  with 
artificial  viscosity  when  the  flow  is  discontinuous.  This  is  important  because  some  fluids 
that  exhibit  spurt  have  negligible  Newtonian  viscosity.  (For  these  fluids,  the  relation 
between  steady  stress  and  strain  rate  resembles  Fig.  2  because  of  multiple  relaxation  time 
scales  [10].) 

The  same  experiment  was  performed  for  system  (C)  with  a  small,  but  nonzero,  New¬ 
tonian  viscosity  coefficient  e.  Figure  8  shows  the  results  for  e  =  0.01,  as  calculated  using 
the  Lax-Wendroff  method  with  Tyler  artificial  viscosity.  What  results  is  a  different  phe¬ 
nomenon,  in  which  the  shorter  relaxation  response  of  the  fluid  (here  modeled  by  Newto¬ 
nian  viscosity)  arrests  the  acceleration  in  a  layer  near  the  wall.  Now  the  slip  layer  is  much 
thicker,  with  its  outer  boundary  corresponding  to  a  discontinuity  in  the  strain  rate  vx. 
The  solution  approaches  a  steady  state  in  which  vx  is  discontinuous  but  the  total  stress 
T  =  a(vx )  +  evx  is  continuous.  The  steady  state  has  the  same  layer  thickness  as  predicted 
analytically,  but  the  centerline  velocity  is  20%  too  high;  this  is  because  the  centerline  ve¬ 
locity  is  extremely  sensitive  to  the  slope  of  the  velocity  profile  in  the  slip  layer,  which  is 
affected  by  the  artificial  viscosity  in  the  numerical  method.  The  layer  formation  is  crucial 
to  our  interpretation  of  the  spurt  phenomenon. 

More  extensive  experiments  were  performed  using  the  solid  mechanics  algorithm.  For 
example,  the  calculation  of  Fig.  8  was  repeated  using  this  method  and  a  graded  mesh  of 
160  elements;  the  same  layer  thickness  as  shown  in  Fig.  8  was  obtained,  and  the  centerline 
velocity  of  the  long-time  solution  differed  from  the  analytic  prediction  by  about  only  1%. 
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Fig.  7:  Onset  of  slip  for  a  fluid  without  Newtonian  viscosity. 

5b.  Spurt  Phenomenon 

We  used  the  solid  mechanics  approach  to  simulate  the  experiments  of  Vinogradov, 
et  a  1.  with  polyisoprene  [18].  In  these  experiments,  the  sample  fluids  are  labeled  PI-1 
through  PI-8,  ordered  by  increasing  molecular  weight  M;  the  parameters  entering  the 
mathematical  model  were  chosen  to  correspond  to  these  samples.  In  the  system  (JS),  a 
measures  the  relative  importance  of  inertial  and  elastic  effects,  and  e  reflects  the  presence 
of  Newtonian  viscosity.  For  the  experimental  fluids,  however,  the  Newtonian  viscosity  is 
negligible  because  there  is  no  solvent;  instead  the  fluids  exhibit  a  second  relaxation  time 
that  is  much  shorter  than  the  dominant  relaxation  time  A-1.  Nevertheless,  the  effects  of 
a  short  secondary  relaxation  time  are  correctly  modeled  by  the  Newtonian  viscosity  term 
provided  that  e  is  interpreted  as  the  ratio  of  the  relaxation  times  [10]. 

The  following  features  of  the  experimental  fluid  samples  were  used  to  determine  the 
physical  constants: 

(i)  The  elastic  modulus  p  is  independent  of  the  molecular  weight  M. 

(ii)  The  contribution  to  the  zero  shear  viscosity  from  the  dominant  relaxation  time,  p/A, 
varies  over  two  orders  of  magnitude  because  of  the  sensitive  dependence  of  A-1  on  M. 

(iii)  There  is  a  critical  molecular  weight  below  which  the  material  does  not  spurt.  (Samples 
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Fig.  8:  Onset  of  spurt  for  a  fluid  with  Newtonian  viscosity. 

PI-1  and  PI-2  do  not  exhibit  spurt.) 

(iv)  For  samples  PI-3-8,  the  critical  stress  for  onset  of  spurt  is  independent  of  M. 

These  observations,  together  with  the  presumption  that  the  secondary  relaxation  time 
is  independent  of  M,  lead  to  values  for  a  and  e  that  decrease  with  M.  These  values 
are  obtained  readily  from  our  definitions  in  §2  and  the  dimensional  information  given  in 
Ref.  [8].  The  computational  results  are  shown  in  Figs.  9-11. 

Figure  9  shows  the  evolution  of  the  spurt  process  in  time;  the  centerline  velocity  is 
plotted  vs.  time  for  sample  PI-7  with  /  =  1.2.  The  spatial  mesh,  with  a  total  of  640 
elements,  was  graded  to  have  smaller  elements  near  the  wall.  All  simulations  were  carried 
out  using  zero  initial  data.  Not  visible  in  Fig.  9  is  a  “start-up”  transient,  occurring  on  a  time 
scale  of  order  e,  during  which  the  centerline  velocity  spikes  to  a  large  value  corresponding 
to  a  Newtonian  fluid  with  viscosity  e.  Following  this  is  a  “latency”  period  when  the 
velocity  and  the  shear  stress  resemble  a  steady  solution  without  spurt.  This  period  ends  at 
t  =  2.36  (i.e.,  376  s,  dimensionally),  when  a  spurt  solution  begins  to  develop.  Notice  that 
these  three  phases  of  the  dynamics  correspond  to  the  results  of  the  phase-plane  analysis 
of  §3c.  Moreover,  the  latency  time  observed  in  the  full  dynamical  simulations  correlates 
precisely  with  the  value  t  =  2.36  found  by  integrating  the  ordinary  differential  equations; 
it  also  agrees  with  the  predictions  of  Eq.  (3.8),  which  gives  t  =  2.30  (see  Ref.  [10]). 
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Fig.  9:  Centerline  velocity  vs.  time. 

The  calculation  shown  in  Fig.  9  has  not  been  run  long  enough  to  achieve  steady  state. 
Essentially  steady  flow  is  attained  after  about  five  more  time  units;  thus  we  predict  that 
the  whole  dynamic  process  for  this  experimental  sample  takes  about  forty  minutes.  We 
have  run  complete  simulations  for  each  of  the  eight  samples.  Figure  10  shows  the  results,  as 
compared  to  the  data  reported  in  Ref.  [18].  In  these  graphs,  effective  shear  rate  5  is  plotted 
vs.  effective  wall  shear  stress  rwaii  =  fh/ 2;  here  the  effective  shear  rate  is  S  =  4Q/ttR3  for 
a  capillary  of  radius  R  and  S  =  6Q/(wh2)  for  a  slit  die  of  width  w,  Q  being  the  volumetric 
flow  rate  [8], 

Figure  11  shows  the  result  of  simulating  a  loading  sequence  in  which  the  pressure 
gradient  /  is  increased  in  small  steps,  allowing  sufficient  time  between  steps  to  achieve 
steady  flow  [8].  The  loading  sequence  is  followed  by  a  similar  unloading  sequence,  in 
which  the  driving  gradient  is  decreased  in  steps.  The  initial  step  used  zero  initial  data, 
and  succeeding  steps  used  the  results  of  the  previous  step  as  initial  data.  The  resulting 
hysteresis  loop  resembles  the  “shape  memory”  observed  by  Hunter  and  Slemrod  [4]  in  their 
model  system  (3.5).  The  width  of  the  hysteresis  loop  at  the  bottom  can  be  related  directly 
to  the  molecular  weight  of  the  sample  [8]  (see  §6). 
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Fig.  10:  Volumetric  flow  rate  vs.  effective  shear  stress:  (a)  experi¬ 
ment  [18];  (b)  numerical  calculation  [8].  Notice  that  the  horizontal 
scale  of  this  panel  matches  that  of  panel  (a),  but  the  vertical  scale 
does  not. 

5c.  Performance  of  the  Numerical  Methods 

Of  the  numerical  methods  we  have  proposed,  the  random  choice  method  is  the  only 
one  that  applies  to  the  interesting  case  in  which  e  =  0;  the  physical  ramifications  of  our 
results  in  this  case  remain  to  be  elucidated.  When  e  ^  0,  the  solid  mechanics  approach 
yields  the  most  fully  developed  and  flexible  method.  The  parabolic  method  may  have 
advantages,  but  we  have  not  attempted  to  refine  it;  our  implementation  is  as  simple  as 
possible,  relying  heavily  on  packaged  codes,  so  as  to  confirm  the  validity  of  the  other  two 
methods. 

Computations  in  the  experimental  regime  yield  interesting  insights  into  the  behavior 
of  our  numerical  methods;  a  case  in  point  is  the  phenomenon  of  latency.  To  accurately 
reproduce  the  latency  time  predicted  by  phase-plane  analysis,  At  must  be  relatively  small 
compared  to  e  during  the  transient  “Newtonian”  phase,  which  takes  place  on  a  time-scale 
of  order  e.  When  the  time  step  is  larger,  the  solid  mechanics  algorithm  remains  stable,  but 
accuracy  is  sacrificed;  thus  the  normal  stress  difference  Z  at  the  end  of  the  Newtonian  phase 
is  not  calculated  precisely.  Indeed,  if  equal  time  steps  are  used  throughout  the  simulation, 
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the  latency  time  is  computed  to  be  t  —  2.14.  That  this  value  compares  favorably  with 
t  ~  2.36  reflects  the  robustness  of  the  algorithm.  By  contrast,  there  is  no  semi-implicit 
treatment  in  our  version  of  the  parabolic  method,  which  requires  a  time  step  smaller  than 
does  than  the  solid  mechanics  approach. 

We  have  been  careful  to  test  the  validity  of  the  numerical  results  we  report  here.  One 
of  the  questions  we  sought  to  resolve  involves  the  oscillations  evident  in  Fig.  9  during  the 
spurt  process.  In  Ref.  [8],  results  were  reported  on  meshes  much  cruder  than  the  one  used 
to  compute  the  results  of  Fig.  9;  the  oscillations  were  larger  in  amplitude,  and  they  did  not 
diminish  with  refinement  of  time  step.  Figure  9  shows  that  these  oscillations  diminish  with 
refinement  of  the  spatial  grid  size,  suggesting  that  they  are  induced  by  spatial  discretization 
error.  This  conclusion  is  reinforced  by  inspection  of  Fig.  9:  the  larger  oscillations  occur 
somewhat  after  the  onset  of  spurt,  when  the  layer  boundary  has  moved  out  of  the  region 
of  mesh  refinement  near  the  wall.  Eventually,  however,  these  oscillations  are  damped,  just 
as  they  are  in  calculations  with  cruder  meshes.  Further  mesh  refinement  studies  indicate 
that  although  crude  spatial  resolution  can  lead  to  spurious  oscillations  in  spurt  dynamics, 
the  solution  maintains  an  accurate  mean  value  and  approaches  the  correct  steady  state. 
These  results  were  reproduced  using  the  parabolic  formulation  with  a  mesh  refined  to  3072 
equal-sized  cells.  To  obtain  accurate  estimates  of  latency  time  with  this  method,  the  time 
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step  must  be  about  half  of  that  required  by  the  solid  mechanics  method;  if  this  condition 
is  violated,  the  spurt  occurs  prematurely.  When  the  time  step  is  so  refined,  the  results 
obtained  with  both  methods  agree  to  graphical  accuracy. 

Our  numerical  approaches  have  successfully  computed  fully  time-dependent  flows  for 
a  Johnson-Segalman  fluid  at  a  high  Deborah  number.  Thus  they  avoid  the  “high  Weis- 
senberg/Deborah  number  problem,”  at  least  in  one-dimensional  flows.  Currently,  we  are 
investigating  generalizations  of  our  approaches  to  multi-dimensional  problems  of  physical 
interest. 

6.  Rheological  Experiments 

In  Ref.  [8],  several  possible  experiments  are  suggested  that  could  verify  the  inter¬ 
pretation  of  spurt  put  forward  here.  The  most  important  experiment  suggested  is  the 
verification  of  the  molecular-weight  dependence  of  the  widest  point  of  the  hysteresis  loop 
of  Fig.  11. 

The  shape  of  the  hysteresis  loop  is  a  key  feature  predicted  by  our  analysis  and  com¬ 
putations:  the  loop  always  opens  from  the  point  at  which  unloading  begins,  and  there 
is  a  discontinuity  in  slope  during  unloading.  This  is  because  the  solutions  proceed  from 
“top-jumping”  in  Fig.  2,  through  intermediate  convexifications  of  the  curve,  to  “bottom¬ 
jumping,”  causing  the  discontinuity  in  slope.  (This  behavior  is  in  distinct  contrast  to  the 
interpretation  of  Ref.  [12],  where  bottom- jumping  is  always  the  rule  for  steady  spurt  solu¬ 
tions,  and  portions  of  the  loading  path  are  retraced  during  unloading.)  The  critical  wall 
stress  for  onset  of  spurt  under  loading  is  Tcr;t  =  1/2  -f  0(e),  whereas  spurt  does  not  cease 
when  unloading  until  the  wall  stress  drops  below  2y/e  [1  +  0(e)].  The  sensitive  dependence 
of  e  on  the  molecular  weight  leads  to  observable  variation  in  the  hysteresis  loop. 

The  analysis  and  computations  presented  in  this  paper  allow  us  to  say  more  about 
experimental  signatures:  (i)  before  dramatic  growth  in  throughput  occurs,  very  slow  flow 
with  little  throughput  persists  during  a  latency  period  that  can  last  several  minutes;  (ii)  la¬ 
tency  occurs  only  wh  u  e  is  sufficiently  small  (i.e.,  less  than  1/8  for  Johnson-Segalman;  the 
precise  number  for  other  models  may  be  different);  (iii)  latency  occurs  only  if  the  applied 
pressure  gradient  is  not  too  large  (/cr it  <  /  <  2  for  Johnson-Segalman);  (iv)  the  latency 
time  scales  with  A-1  when  the  pressure  gradient  is  fixed,  and  is  determined  approximately 
by  integrating  Eq.  (3.8). 

7.  Conclusions 

Well-posed  dynamical  problems  for  fluids  with  non-monotone  constitutive  relations 
need  not  be  unphysical.  In  fact,  the  Johnson-Segalman  model  provides  a  relatively  simple 
example  that  accurately  describes  spurt.  Our  analytical  and  computational  results  predict 
several  phenomena  related  to  spurt,  which  should  be  observable  in  rheological  experiments. 
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